Inversion#

In what follows, we will implement inversion methods for sampling discrete and continuous distributions.

Discrete distributions#

Example 2.1. We think of our discrete distribution very practically as a vector defined on a set of states. Let us define a probability mass function on the set \(\mathsf{S} = \{1, 2, 3, 4, 5\}\) as follows:

\[\begin{align*} \mathsf{w} = \begin{bmatrix} 0.2 & 0.3 & 0.2 & 0.1 & 0.2 \end{bmatrix} \end{align*}\]

Let us plot the PMF and CDF of this distribution.

import numpy as np
import matplotlib.pyplot as plt

w = np.array([0.2, 0.3, 0.2, 0.1, 0.2])
s = np.array([1, 2, 3, 4, 5])

def discrete_cdf(w):
    return np.cumsum(w)

cw = discrete_cdf(w)

def plot_discrete_cdf(w, cw):
    fig, ax = plt.subplots(1, 2, figsize=(20, 5))
    ax[0].stem(s, w)
    ax[1].plot(s, cw, 'o-', drawstyle='steps-post')
    plt.show()

plot_discrete_cdf(w, cw)
../_images/1bdf8d16d42a9725d23f10a3c863358f24e64aac059004a3f9b90710fbd8a883.png

Next, we can implement the sampling method using the inverse. Let us implement the method first, then look into some animations.

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(4)

w = np.array([0.2, 0.3, 0.2, 0.1, 0.2])
s = np.array([1, 2, 3, 4, 5])

def discrete_cdf(w):
    return np.cumsum(w)

def sample(u, s, w): # discrete sampler for a uniform random variable u, states s and weights w
    cdf = discrete_cdf(w)
    sample_ind = np.argmax(cdf > u)
    return s[sample_ind]

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
n = 2000
x = []
cw = discrete_cdf(w)
un = []

for n in range(n):
    u = rng.uniform(0, 1)
    un.append(u)
    sample_x = sample(u, s, w)
    x.append(sample_x)

    # the rest is for animation

# plot the pmf using the stem function and change the color of the markers and the line
markerline, stemlines, baseline = ax[0].stem(s, w, markerfmt='o', linefmt='r-')
ax[0].cla()
ax[0].hist(un, bins=10, density=True, color='k', alpha=1)
ax[0].set_title("Histogram of Uniform Random Variables")

ax[1].cla()
ax[1].stem(s, w, markerfmt='o', linefmt='r-')
ax[1].set_xlabel("s")
# plot a histogram centered on states s
ax[1].hist(x, bins=range(7), density=True, color='k', alpha=1, align='mid', width=0.1)
ax[1].set_xlim([0, 6])
ax[1].set_title("Histogram of Discrete Random Variables and PMF")
plt.show()
../_images/263943cefa8330daee56299f489fea2c3f18f69cd30b2ab6cd894c60c4c09823.png

As we can see, with \(n = 2000\), the method samples from the correct distribution. Let us see below this process animated (the code is hidden as it has a lot of diversions from normal code for the sake of animation – but feel free to expand if you are curious!).

Hide code cell source
%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display

rng = np.random.default_rng(4)

w = np.array([0.2, 0.3, 0.2, 0.1, 0.2])
s = np.array([1, 2, 3, 4, 5])

def discrete_cdf(w):
    return np.cumsum(w)

def sample(u, s, w): # discrete sampler for a uniform random variable u, states s and weights w
    cdf = discrete_cdf(w)
    sample_ind = np.argmax(cdf > u)
    return s[sample_ind]

fig, ax = plt.subplots(2, 2, figsize=(10, 8))
n = 200
x = []
cw = discrete_cdf(w)
un = []

def update(i):
    global un, x
    u = rng.uniform(0, 1)
    un.append(u)
    sample_x = sample(u, s, w)
    x.append(sample_x)

    # the rest is for animation

    if i % 1 == 0:
        ax[0, 0].cla()
        # plot the pmf using the stem function and change the color of the markers and the line
        markerline, stemlines, baseline = ax[0, 0].stem(s, w, markerfmt='o', linefmt='r-')
        ax[0, 0].set_title("PMF")
        ax[0, 0].set_xlabel("s")
        # plot u in the y axis of ax[0, 1]
        ax[0, 1].cla()
        ax[0, 1].plot(s, cw, 'ro-', drawstyle='steps-post')
        ax[0, 1].set_title("Cumulative Distribution Function")
        ax[0, 1].set_xlabel("s")
        ax[0, 1].set_xlim([0, 6])
        ax[0, 1].set_ylim([0, 1])
        ax[0, 1].plot(0, u, c='k', marker='o', linestyle='none', markersize=10)
        ax[0, 1].plot(sample_x, u, c='k', marker='o', linestyle='none', markersize=10)
        ax[0, 1].plot([0, sample_x], [u, u], c=[0.8, 0, 0], linestyle='--')
        ax[0, 1].plot(sample_x, 0, c='k', marker='o', linestyle='none', markersize=10)

        ax[1, 0].cla()
        ax[1, 0].hist(un, bins=10, density=True, color='k', alpha=1)
        ax[1, 0].set_title("Histogram of Uniform Random Variables")

        ax[1, 1].cla()
        ax[1, 1].stem(s, w, markerfmt='o', linefmt='r-')
        ax[1, 1].set_xlabel("s")
        # plot a histogram centered on states s
        ax[1, 1].hist(x, bins=range(7), density=True, color='k', alpha=1, align='mid', width=0.1)
        ax[1, 1].set_xlim([0, 6])
        ax[1, 1].set_title("Histogram of Discrete Random Variables and PMF")


ani = FuncAnimation(fig, update, frames=n, repeat=False)
HTML(ani.to_jshtml())

Of course the idea generalises to continuous distributions. Now let us demonstrate Example 2.2 of the lecture notes.

Example 2.2. In the lecture we have derived the sampler for exponential distribution:

\[\begin{align*} p(x) = \text{Exp}(x;\lambda) = \lambda e^{-\lambda x}. \end{align*}\]

We calculate the CDF

\[\begin{align*} F_X(x) =& \int_{0}^x p(x') \mathrm{d} x', \\ =& \lambda \int_0^x e^{-\lambda x'} \mathrm{d} x', \\ =& {\lambda} \left[ -\frac{1}{\lambda} e^{- \lambda x'}\right]_{x' = 0}^x \\ =& 1 - e^{-\lambda x}. \end{align*}\]

Deriving the inverse:

\[\begin{align*} u =& 1 - e^{-\lambda x}\\ \implies x =& -\frac{1}{\lambda} \log (1 - u)\\ \implies F_X^{-1}(u) =& -\lambda^{-1} \log (1 - u). \end{align*}\]

which gives us the sampler:

  • Generate \(u_i \sim \text{Unif}([0, 1])\)

  • \(x_i = -\lambda^{-1} \log (1 - u_i)\).

Let us look at the code for this sampler.

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(10)

def exponential_pdf(x, lam):
    return lam * np.exp(-lam * x)

def exponential_cdf(x, lam):
    return 1 - np.exp(-lam * x)

# an illustration of inverse transform sampling for the exponential distribution
# with parameter lambda = 4
lam = 4
N = 50000

x = np.array([])
un = []
for n in range(N):
    u = rng.uniform(0, 1)
    un.append(u)
    x = np.append(x, -np.log(1 - u) / lam)

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
# plot u on the y axis
xx = np.linspace(0, 2, 100)
axs[0].hist(un, bins=50, density=True, color='k')
axs[0].set_title("Histogram of Uniforms")
axs[1].hist(x, bins=50, density=True, color='k')
axs[1].plot(xx, exponential_pdf(xx, lam), color=[0.8, 0, 0], linewidth=2)
axs[1].set_xlim(0, 2)
axs[1].set_title("Histogram of Exponentials")
plt.show()
../_images/eb385e755efcfda4b58f41b280082d1cedd4f14cb4ddbd8169581b1b19b8becc.png

We can see that the sampler works as expected. Let us animate the process.

Hide code cell source
%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display

rng = np.random.default_rng(4)

def exponential_pdf(x, lam):
    return lam * np.exp(-lam * x)

def exponential_cdf(x, lam):
    return 1 - np.exp(-lam * x)

# an illustration of inverse transform sampling for the exponential distribution
# with parameter lambda = 4
lam = 4
N = 200

fig, axs = plt.subplots(2, 2, figsize=(7, 7), width_ratios=[1, 3], height_ratios=[3, 1])
x = np.array([])
un = []
def update(n):
    global un, x
    u = rng.uniform(0, 1)
    un.append(u)
    x = np.append(x, -np.log(1 - u) / lam)

    # print(x[n])
    if n % 1 == 0:
        # plot u on the y axis
        axs[0, 0].cla()
        axs[0, 1].cla()
        xx = np.linspace(0, 2, 100)
        axs[0, 1].plot(xx, exponential_cdf(xx, lam), color='k', linewidth=2)
        axs[0, 1].plot(0, u, 'o', color=[0.8, 0, 0], markersize=4)
        # plot the projection of u on the pdf
        axs[0, 1].plot(x[-1], u, 'o', color=[0.8, 0, 0], markersize=4)
        # plot a dashed line between the two points
        axs[0, 1].plot([0, x[-1]], [u, u], '--', color=[0.8, 0, 0], linewidth=2)
        axs[0, 1].plot([x[-1], x[-1]], [u, 0], '--', color=[0.8, 0, 0], linewidth=2)
        axs[0, 1].set_xlim(0, 2)
        axs[0, 1].set_ylim(0, 1)
        axs[0, 1].plot(x[-1], 0, 'o', color=[0.8, 0, 0], markersize=4)
        # plot on the axis on the left the histogram of un parallel to the y axis of the axs[0, 1] plot
        axs[0, 0].hist(un, bins=50, density=True, color='k', orientation=u'horizontal')
        axs[0, 0].set_ylim(0, 1)
        axs[0, 0].set_xlim(0, 1)
        axs[0, 0].invert_xaxis()
        axs[1, 1].cla()
        axs[1, 1].hist(x, bins=50, density=True, color='k')
        axs[1, 1].plot(xx, exponential_pdf(xx, lam), color=[0.8, 0, 0], linewidth=2)
        axs[1, 1].set_xlim(0, 2)
        # set horizontal orientation of the histogram on the right
        axs[1, 0].set_visible(False)

ani = FuncAnimation(fig, update, frames=N, repeat=False)
HTML(ani.to_jshtml())

One can see above the animation of uniform sample being sampled and going through the inverse of the CDF which eventually becomes a sample from the exponential density.

Using these examples, you can practice to implement other inversion methods in lecture notes. A good idea is to animate other examples using the code above.